Scaling, Multiscaling, and Nontrivial Exponents in Inelastic Collision Processes 
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We investigate velocity statistics of homogeneous inelastic gases using the Boltzmann equation. 
Employing an approximate uniform collision rate, we obtain analytic results valid in arbitrary di- 
mension. In the freely evolving case, the velocity distribution is characterized by an algebraic large 
velocity tail, P(v,t) ~ v~ a . The exponent a(d, e), a nontrivial root of an integral equation, varies 
continuously with the spatial dimension, d, and the dissipation coefficient, e. Although the velocity 
distribution follows a scaling form, its moments exhibit multiscaling asymptotic behavior. Further- 
more, the velocity autocorrelation function decays algebraically with time, A(t) = (v(0) •v(t)) ~ t~ a , 
with a non-universal dissipation-dependent exponent a = 1/e. In the forced case, the steady state 
Fourier transform is obtained via a cumulant expansion. Even in this case, velocity correlations 
develop and the velocity distribution is non-Maxwellian. 
PACS: 05.20.Dd, 02.50.-r, 47.70.Nd, 45.70.Mg 



I. INTRODUCTION 

Inelastic gases consist of hard sphere particles that 
interact via contact interactions and dissipate kinetic 
energy upon collisions pj. They are used extensively 
to study dynamics of granular materials. Numerically, 
molecular dynamics simulations are quite successful in 
modeling many of the observed collective phenomena 
that include size segregation, phase transitions, shocks, 
clustering, and development of other spatial structures 
[p|-|lO| . In parallel, kinetic theory is utilized to systemat- 
ically derive macroscopic properties from the microscopic 
collision dynamics |Il"|-|l3|] . 

Inelastic gases, a prototype noncquilibrium interact- 
ing particle system, are interesting on their own rights 
p4j-pl[. Recent theoretical and experimental studies 
show that the velocity distributions exhibit anomalous 
large-velocity statistics with exponential, stretched ex- 
ponential, and Gaussian tails |,||^-|25| . Inelastic gases 
involve significant velocity and spatial correlations in 
contrast with traditional molecular gases |p6|-p8|. Ki- 
netic theory assumes that spatial velocity correlations 
are small. While this assumption can be justified for 
strongly driven gases, the situation for freely evolving 
gases is more difficult since velocity correlations can be 
ignored only in the early homogeneous phase p^-pl[| , but 
must be taken into account in the asymptotic clustering 
phase. Clearly, the strong energy dissipation raises chal- 
lenging new questions Q . 

Yet, even more elementary questions remain unan- 
swered. For example, random collision processes effec- 
tively generate thermal, purely Maxwellian, velocity dis- 
tributions when the collisions are elastic. In particu- 
lar, different components of the velocity become uncor- 
rected. In this study, we consider these very same pro- 
cesses but with inelastic collisions. We show that energy 
dissipation fundamentally alters the behavior. The sys- 
tem is intrinsically a noncquilibrium one, and the result- 
ing velocity distributions are far from thermal. 

We consider a collision process where random pairs of 



particles undergo inelastic collisions with a random im- 
pact direction. This process, often called the Maxwell 
model, is described by a Boltzmann equation with a uni- 
form collision rate. In classical kinetic theory of gases, 
the analytically tractable Maxwell model precedes the 
Boltzmann equation | j33[ . Historically, it played an im- 
portant role in the development of kinetic theory [ p4| 36 
and it still remains the subject of current research p7| . 

Very recently, it has been noted that the Maxwell 
model is analytically tractable even for inelastic colli- 
sions [p9fji^|. Interesting behavior emerges in the freely 
evolving case. In one dimension, while moments of the 
velocity distribution exhibit multiscaling |59|, the veloc- 
ity distribution itself still approaches a scaling form with 
an algebraic large velocity tail pQ] . Here, we show an- 
alytically that in arbitrary spatial dimension the veloc- 
ity distribution admits a scaling solution with an alge- 
braic large velocity tail. The corresponding exponent, a 
root of a transcendental equation, depends on the spa- 
tial dimension and the restitution coefficient. Addition- 
ally, we find that the multiscaling behavior extends to 
higher dimensions, and that the velocity autocorrela- 
tion function exhibits aging and nonuniversal asymptotic 
behavior. In general, velocity components develop sig- 
nificant correlations. Such correlations diminish in the 
forced case, although the velocity distribution remains 
non-Maxwellian . 

The rest of this paper is organized as follows. The 
basic Boltzmann equation for the velocity distribution 
and its Fourier transform are presented in Sec. II. In 
Sec. Ill, we investigate the scaling regime, and obtain 
the extremal velocity statistics, moments of the velocity 
scaling function, and velocity correlations. In Sec. IV, we 
illuminate the nonequilibrium dynamics by studying the 
time dependent behavior of the moments and the velocity 
autocorrelation function. In Sec. V, we consider nonequi- 
librium steady states in the driven case, and obtain the 
steady state distribution as a cumulant expansion. A 
few generalizations are briefly mentioned in Sec. VI, and 
conclusions are given in Sec. VII. 
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II. THE MAXWELL MODEL 

We study a homogeneous system of identical inelastic 
spherical particles. The mass and the cross-section are 
set to unity without loss of generality. Particles interact 
via binary collisions that lead to exchange of momentum 
along the impact direction. The post-collision velocities 
Vi ; 2 are given by a linear combination of the pre-collision 
velocities u 1-2 , 



vi, 2 = ui )2 T (1 - e) (g • n) n. 



(1) 



Here g = ui — 112 is the relative velocity and n the unit 
vector connecting the particles' centers. In each collision, 
the normal component of the relative velocity is reduced 
by the restitution coefficient r = 1 — 2e. The energy 
dissipation equals AE = — e(l — e)(g • n) 2 , so for e = 
collisions are elastic, while for e = 1/2 collisions are per- 
fectly inelastic with maximal energy dissipation. Since 
the collision rule ([!]) is Galilean invariant, the average 
velocity can be set to zero without loss of generality. 

We investigate the "Maxwell model" where the colli- 
sion rate in the Boltzmann equation equals the typical 
velocity, rather than the actual relative velocity |p5| , ^6| . 
This kinetic theory describes a stochastic process where 
randomly chosen pairs of particles undergo inelastic col- 
lisions according to ([!]) with a randomly chosen impact 
direction n. In such a process, no spatial correlations de- 
velop, and the normalized velocity distribution function, 
P(v,t), obeys 



dP(y,t) 



. )t g J dnj duiP(ui,t) / du 2 P(u 2 ,t) (2) 
x {5 [v - ui + (1 - e)(g • n)n] - S(y - u x )}. 



The overall collision rate equals g = VT where T is the 
granular temperature, or the average velocity fluctuation 
per degree of freedom, T = ^ J dw 2 P(v,t) with v = |v|. 
The restriction g • n > on the angular integration range 
in Eq. can be tacitly ignored, because the integrand 
obeys the reflection symmetry n ► — n. This angular 
integration should be normalized, J dn = 1. 

We study primarily the freely evolving case where in 
the absence of energy input the system "cools" indefi- 
nitely. From the Boltzmann equation (Q), the tempera- 
ture rate equation is 



—T = —AT 3 / 2 , with 
dt 



The constant A = 2e(l — e) J dnnf, is obtained using the 
identity n\ 



Given the convolution structure of the Boltzmann 
equation we introduce the Fourier transform [[3G| of 
the velocity distribution function, 



F(k,t) = J dve ikv P(v,i). 



(5) 



We conveniently reset the collision rate to unity by modi- 
fying the time variable. The collision counter r is defined 
via the transformation = ^7=^- Specifically, 



•*/** 



(6) 



equals the average number of collisions experienced by a 
particle. Applying the Fourier transform to Eq. (^|) and 
integrating over the velocities gives 



P(k,r) 



dnP[k-q,r]F[q,r], (7) 



with q = (1 — e)(k • n) n. This equation reflects the mo- 
mentum transfer occurring during collisions. 

We restrict our attention to isotropic situations, and 
write the Fourier transform P(k, t) = F(z, t) in terms of 
the variable z — k 2 . To perform the angular integration, 
it proves useful to employ spherical coordinates with the 
polar axis parallel to k, so that k • n = cos#. The 9- 
dependent factor of the measure dn is proportional to 
(sin 9) d ~ 2 d9. In terms of the variable /1 = cos 2 9 one has 
dn = T>fi, with 



2' 2 



V/i = i-i 2 (1 - n) 2 dfi 



(8) 



where B(a, b) is the beta function. This integration mea- 
sure is properly normalized, J T>fi = 1. Hereinafter, we 
denote angular integration with brackets 



(/> = / Z> M /(/*)• 



(9) 



The governing equation (Q) for the Fourier transform can 
now be rewritten in the convenient from 



d_ 



F(z, t) + F(z, t) = ( F(£z, t) F( V z, t) ), (10) 



2 - 1 that vields frinn 2 - Md with the shorthand notations £ 
d - 1 tnat yields J an ni - i/cl. _ 2 v™,™ 



1-(1 



)n and 



Solving Eq. (||) we find that the temperature decays ac- 
cording to Haff 's cooling law [Fbil 



T(t)=T (l + t/Q- 



(4) 



with the time scale £* = dj [e(l — e)V7o] set by the ini- 
tial temperature, Tq. 



i] = (1 — e) 2 [i. Hence, the Fourier equation is both non- 
linear and non-local. Interestingly, while it is difficult to 
integrate this equation with respect to time, most of the 
physically relevant features of the velocity distributions 
including large velocity statistics and the time depen- 
dent behavior of the moments can be found analytically, 
as will be shown below. 
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III. SCALING SOLUTIONS 

Numerical simulations in two-dimensions suggest that 
the velocity distribution approaches the scaling form jit]] 



P(v,i) 



T d/2' 



The scaling form of the Fourier transform reads 



F(k,i) = $(x), with 



k 2 T. 



(11) 



(12) 



In the k — > limit, the Fourier transform behaves as 
F(k,t) = 1 - \k 2 T. This implies that the first two 
terms in the Taylor expansion of the corresponding scal- 
ing function are universal, $(x) = 1 — \x. Substitut- 
ing the above scaling form into the governing equation 
( |l0| ) and using the temperature cooling rate -§^T = —XT 
yields the governing equation for the scaling function 



Ax $'(x) + $(x) = (<I>(£x) $(r?x)). 



(13) 



One can check that the velocity distribution is purely 
Maxwellian <E>(x) = e~ x l 2 in the elastic case [Q. Indeed, 
A = and £ + r\ = 1 in this case. A stochastic process of 
elastic collisions effectively randomizes the velocities and 
leads to a thermal distribution. 



A. Algebraic tails 

It is instructive to consider first the one-dimensional 
case. Here, integration over fi is immediate as this 
variable equals unity, and the scaling function satisfies 
-Ax<j>'(x) + <I>(x) = $ [e 2 x] $ [(1 - e) 2 x\. Remarkably, 
this non-local non-linear differential equation admits a 
very simple solution |if]] 



$(x) = (l + y/x) 



(14) 



Performing the inverse Fourier transform gives the veloc- 
ity distribution as a squared Lorentzian 



V(w) 



2 
n 



1 



w 



(15) 



The scaling solution (|15|) is universal as it is indepen- 
dent of the dissipation coefficient e. Its key feature is the 
algebraic tail, "P(u>) ~ w~ 4 as \w\ — > oo. 

In general dimension d, the large velocity behavior 
of the velocity distribution can be determined from the 
small wave number behavior of its Fourier transform. For 
example, the small-x expansion of the one-dimensional 
solution (|l4|) contains both regular and singular terms: 
$(x) = 1 — | x + | x 3 / 2 H , and the dominant singu- 
lar x 3 / 2 term reflects the w~ 4 tail of V(w). In general, 
an algebraic tail of the velocity distribution (|l] 



V(w) 



as w 



oo, 



(16) 



indicates the existence of a singular component in the 
Fourier transform, 



0. 



(17) 



The inverse is also correct. This can be seen by recasting 
the Fourier transform <E>(x) oc dww d ~ 1 V(w) e lw ^ 
into a Laplace transform I(s) oc J °° dw w d ~ 1 V(w) e~ ws 
by writing x = — s 2 . The small- s expansion of I(s) 
contains regular and singular components. For exam- 
ple, when a < d, the integral I(s) diverges as s — * 
and integration over large-w yields the dominant contri- 



buti> 



'sing 



00 



When d < a < d + 1, 1(0) is 



finite, but the next term is the above singular term, so 



1(a) = 1(0) 



■. In general, the singular con- 
d , thereby leading to Eq. (fl7l). 
The exponent a can be now obtained by inserting 



tribution is ising(s) 



$(x) = $ rog (x) + < &sing(x) into Eq. ( |13j ) and balancing 
the dominant singular terms. We find that a is a root of 
the integral equation 



1 — A - 



{a-d)/2 



V 



{a 



-d)/2\ 



(18) 



This relation, originally derived in Refs. [^ljj42j], can 
be recast as an eigenvalue problem. Indeed, defining 
we can re-write Eq. as A M = /1A1 



1 



-d 



with /i = ^j^. Note that A = Ai, so there is an obvious 
solution // — 1, or a = d + 2. In this case, the singular 
term simply coincides with the dominant regular term, 
x {v-d)/2 _ x ^ jj encej solution is trivial and in the 
following we shall seek a solution with a > d + 2. 

The integral equation (|l^) can also be rewritten in 
terms of special functions. The first integral on the right- 
hand side of (|l8|) can be expressed in terms of the hyper- 
geometric function 2-F1 (a, b; c; z) p|| and the second as a 
ratio of beta functions: 



1 - 
2F1 



d- 



' d 
aid 
~ 1 2 ; 2 ; 



(19) 



_ d r(^|±i)r( 



r(f)r(|) 



We conclude that the exponent a = a(d, e) depends in a 
nontrivial fashion on the spatial dimension d as well as 
the dissipation coefficient e. 

First, let us investigate the dependence on the dissi- 
pation coefficient by considering the quasi-elastic limit 
e — > 0. In the elastic case, the Maxwellian distribution, 
$(x) = e~ x / 2 , implies a diverging exponent a — > 00 as 
e — > 0. Therefore, the right-hand side of Eq. ([l9]) van- 
ishes in the quasi-elastic limit and to leading order 



d 

a ~ — . 

e 



(20) 



Clearly, the quasi-elastic limit is singular. Dissipation, 
even if minute, seriously changes the nature of the sys- 
tem |||5 31 1 . Further corrections can be obtained via a 



3 



systematic perturbation expansion in e. We merely quote 
the two leading corrections in the physically relevant di- 
mensions 



a(2,e) 
<r(3,e) 



2 (e~ 2 + 1) 4tt - (e- 2 + if 




0(e 1/2 ), 



6-e 



■Oie 1 / 2 ) 



Next, we discuss the dependence on the dimension. 
First, one can verify that a = 4 when d = 1 by utiliz- 
ing the identity 2-Fi(a, b; b; z) = (1 — z)~ a . In the in- 
finite dimension limit, the second integral (j^* 7- ^/ 2 ) in 
Eq. ( |l8| ) is negligible as it vanishes exponentially with 
the dimension. To evaluate the second integral we take 
the limits d — > oo and \i — > with z = /i<i/2 fixed. 
The integration measure (||) is transformed according to 
— > (7rz) _1 / 2 e -2 dz, and the basic Eq. ( |l8| ) becomes 
e(l - e )„ = J oo dz(^z)- 1 /2 e -[i+(i- e 2 H- with thc 



1 

shorthand notation u 
tion yields 1 — e(l — e)u 



— 1. Performing the integra- 
[l + (l-e 2 )u]- 1/2 . This cubic 
equation has the aforementioned trivial solution u = 
and two non-trivial solutions. Choosing the physically 
relevant u, we obtain that as d — > oo 



a 
d 



1 



2 C c 



,1/2 



(1 



,1/2 



6(1 -6 2 ) 



(21) 




FIG. 1. The exact exponent a, obtained from Eq. (|l9|), ver- 
sus the dissipation parameter e. The exponent was scaled by 
the dimension d. Shown also is the limiting large dimension 
expression (pH). 

In general, a oc d, and therefore, the algebraic decay 
becomes sharper as the dimension increases. The ex- 
ponent <j(d, e) increases monotonically with increasing d, 
and additionally, it increases monotonically with decreas- 
ing e (see Fig. 1). Both features are intuitive as they 
mirror the monotonic dependence of the energy dissipa- 
tion rate A = 2e(l — e)/d on d and e. Hence, the com- 



pletely inelastic case provides a lower bound for the ex- 
ponent, a(d, e) > a(d, e = 1/2) with <r(d, 1/2) = 6.28753, 
8.32937, for d = 2, 3, respectively. Numerical simulation 
results are consistent with the former value Thc 
algebraic tails are characterized by unusually large ex- 
ponents which may be difficult to measure accurately in 
practice; for instance, typical granular particles are char- 
acterized by the dissipation coefficient e w 0.1 yielding 
a w 30 in three dimensions. Figure 1 also shows that 
the quantity a/d weakly depends upon the dimension, 
and the large-d limit ( ^l| ) provides a good approximation 
even at moderate dimensions. 



B. Divergence of the moments 

The algebraic tail of the velocity distribution implies 
that sufficiently small moments of the scaling function 
<f>(x) are finite, while moments larger than some index 
diverge. In the scaling regime, moments of the velocity 
distribution can be calculated by expanding the Fourier 
transform in powers of x, 



E 

n>0 



(22) 



The coefficients <j> n yield thc leading asymptotic behav- 
ior of the velocity moments, Mk(t) — J dvv k P(v,t), via 
the relation (2n)!T™0„ ~ (/x")M2 n . Inserting the mo- 
ment expansion into the governing equation ([l3]) yields 
the closed hierarchy of equations 



(A n - nXl)<j)n = 



n-1 

£■ 

m— 1 



^m.n—m^mH^n—mi 



(23) 



with A„ = (l-e i -V n ) and X m ,i = (fV)- Th e first few 
coefficients are written explicitly in Appendix A. Starting 
with cf>o = 1 and <f>\ — 1/2, further coefficients are deter- 
mined recursively from (p3|). In the elastic case (e = 0), 
one has <fi n = (n!2")~ 1 , consistent with <&(x) = e~ x / 2 . 
For general e, the first two terms are 



11-3 



d+2 



1 - 3 



1 1 



(24) 



d+2 

■ 1-e 2 
' d+2 



o (l-e)(l+3e) , on e(l-e)(l-e 2 ) 
d+2 ~T OU (d+2)(d+4) 



1-3 



(i-h) 2 

d+2 



10 



;(l-6)(3+62) 
(d+2)(d+i) 



The behavior is determined by two parameters: d and e. 
Fixing e, we see that a given moment <f> n is finite only if 
the dimension is sufficiently large, d > d n (e). In partic- 
ular, 4> n is finite only if the left-hand side of Eq. ( p3| ) is 
positive, A„ — nAi > 0. This condition is satisfied only if 
the dimension is sufficiently large d > d n , with d n being 
the spatial dimension at which A„ = nAi. For example, 
4>2 > when d > c?2, and (f>3 is finite only when d > d% 
with thc following crossover dimensions 



4 



d 3 



f 3e 2 , 

(e 2 + 2e - 1 



(25) 



60e + 186e 2 - 4e 3 + 49e 4 



Conversely, for a fixed dimension, a given moment is fi- 
nite only if the dissipation is sufficiently small. For ex- 
ample, fa is positive only when e < 0.302074,0.427438 
at d — 2,3. Such values, obtained by solving polynomial 
equations yield integer values of the large velocity decay 
exponent a(2, 0.302074) = 8 and ct(3, 0.427438) = 9, in 
accord with direct numerical solution of Eq. fll9| ) . 



A. Multiscaling of the moments 

While moments of the scaling function diverge, the ac- 
tual moments must remain finite at all times, particu- 
larly at the scaling regime. Therefore, the above moment 
analysis suggests that knowledge of the leading asymp- 
totic behavior is not sufficient to characterize the time 
dependent behavior of sufficiently large moments. 

The time evolution of the moments can be studied us- 
ing the expansion 



C. Velocity correlations 



F(z,r)=J2Ur) (-z)> 



(28) 



Maxwell's seminal derivation of the Maxwellian distri- 
bution (see Ref. j26]| , p. 36) relies on two basic assump- 
tions: (1) Isotropy of the velocity distribution, and (2) 
Absence of correlations between the velocity components. 
The latter assumption is directly probed using the follow- 
ing correlation measure 



Q 



y' 



(26) 



A non-vanishing Q indicates that velocity correlations 
do exist, and the larger Q the larger the correlation. 
In the freely evolving case, this quantity easily fol- 
lows from the small- X behavior of the scaling function 
<$>(x). By definition, (u 2 ) = (u 2 ) = T and furthermore, 



( v 2 v 2 ) - -^jr^F 

{VxV v> k=o 
has Q = 4$"(0) -1 = 8 



= 4T 2 $"(0). Consequently, one 
6 2 -l. Using Eq. (§3) we find 



6e 2 



d-(l + 3e 2 )' 



(27) 



when d > d2 = 1 + 3e 2 , and Q = 00 otherwise. While the 
quantity Q is physical when d > 2, it is sensible to use an- 
alytic continuation to reveal the underlying divergence. 
Velocity correlations vanish for elastic gases. Interest- 
ingly, inelasticity introduces strong velocity correlations, 
and the larger e the larger the correlations as Q increases 
monotonically with increasing e. The perfectly inelastic 
case (e = 1/2) again provides a bound: Q < Q max = 6, 
6/5 for d = 2, 3, respectively. This behavior is somewhat 
intuitive as the unisotropic collision rule (jj]) discriminates 
the velocity component normal to the impact direction. 



IV. NONEQUILIBRIUM DYNAMICS 

Thus far, we focused on the leading asymptotic behav- 
ior of the velocity distribution. The diverging moments 
and the dissipative nature of this system suggest that the 
time dependence may exhibit rich behavior. Thus, we 
study relaxation of velocity characteristics such as the 
moments and the autocorrelation function. 



The actual moments are related to the coefficients via 
(2n)\f n = ([i n )M 2n . Substituting the expansion (pl|) into 
(M) yields the evolution equations 



7 fn "f" ^-nfn 

dr 



^ ^m,n—mfmfn 



(29) 



We demonstrate multiscaling asymptotic behavior by 
evaluating the second, fourth, and sixth moments. The 
second moment is obtained from -rpfi + Xxfx = with 
Ai = A = 2e(l — e)/d. Hence, we recover Haff's 
law /i(r) = /i(0)e- Al -, or fx(t) = fx(0) (1 + t/uy 2 . 
Asymptotically, the second moment of the velocity dis- 
tribution has the universal behavior, M2 ~ t~ 2 . The 
next coefficient f 2 satisfies 



d f 



A2/2 — M,ifi ■ 



(30) 



Solving Eq. ( |30| ) we find that /2(t) is a linear combination 
of two exponentials, e~ A2T and e~ 2Air , whose decay coef- 
ficients are equal A2 = 2Ai at the crossover dimension d 2 . 
Integrating the rate equation (^) and translating back 
to the physical time t, we obtain 



f 2 {t) = C x (1 + t/U) 4 + C 2 (! + £/**) 



-2a 2 



(31) 



ioxd^d 2 . Here, a n = X n /X u d = A M / 1 2 (0)/(A 2 -2Ai), 
and C2 = /2(0) — C\. When d = d 2 one finds 



/ 2 (*) = [Cx In (1 + t/U) + C 2 ] (1 + t/U)' 



(32) 



with Cx = Xx,ifx(°) and C2 = / 2 (0). Thus for d > d 2 , 
the fourth moment exhibits ordinary scaling, f 2 ~ £~ 4 , 
or M A ~ Mf ~ T 
apparent as f 2 ~ 



When d < d 2 , multiscaling becomes 



t~ 2a2 and therefore the ratio M4/M 2 
diverges asymptotically. This is consistent with the di- 
vergence of the fourth moment of the scaling function 
$(x) that occurs at the same crossover dimension d 2 . A 
logarithmic correction occurs at this dimension. In sum- 
mary, we find the following leading asymptotic behavior 
of the fourth moment 



5 



M 4 (t) 



(33) 



t- 4 d>d 2 , 
i~ 4 hit d — d 2 , 
I t- 2a2 d < d 2 . 

A similar calculation can be carried for the sixth mo- 
ment. The solution of £pf 3 + A3/3 = (Ai, 2 + A 2 ,i)/i/2, 
with fx and f 2 given above involves three exponentials: 



-A 3 r 



, e _(Al+A2)r , and e" 



3Ajr 



Asymptotically, the first 



exponential dominates when d < d$, and consequently, 
M 6 ~ t" 2 " 3 ; otherwise, the third exponential dominates 
and thence ordinary scaling occurs, Mq ~ i~ 6 . Generally, 
the leading asymptotic behavior of the 2n-th moment is 
characterized by two different regimes 

(t -2n 



M 2n {t) 



t- 2n lni 

t -2a n 



d > d n 
d = d n 
d < d n 



(34) 



Further logarithmic corrections affecting sub-dominant 
terms occur at the crossover dimensions d 2 , . . . , d n -\. 

The dependence of d n (t) on the dissipation coefficient 
is shown in Figure 2. In the physical dimensions d = 2, 3, 
the fourth moment exhibits ordinary scaling behavior. 
The sixth order moment exhibits multiscaling if the dissi- 
pation coefficient is large enough: e > 0.302074, 0.427438 
for d = 2, 3, respectively. In the large n limit, A„ — I 
so from A n = n\\ we find d n —> 2e(I — e) n. Thus, re- 
gardless of the dissipation parameters and the dimension, 
sufficiently large moments exhibit multiscaling: 



M 2n oc Af 9 c 



A n /Ai. 



(35) 



Interestingly, the multiscaling exponents saturate asymp- 
totically, a n — > d/[2e(l — e)] as n — > 00. Of course, if 
the dimension increases or the dissipation parameter de- 
creases, the order of the lowest moment exhibiting mul- 
tiscaling increases, and in practice, it may be difficult 
to observe deviations from ordinary scaling. For exam- 
ple, at d = 3 and e = 0.1, multiscaling occurs only for 
moments whose index exceeds 30! 




FIG. 2. The crossover dimensions d n (e) of Eq. (|2E]) versus 
the dissipation coefficient for n = 2, 3. 



B. Non-universal velocity autocorrelations 

The autocorrelation function quantifies memory in the 
velocity of a tagged particle j26) . The velocity autocor- 
relation, A(t w ,t), is defined via 



A(t w ,t) = v(t w )-v(t) 



(36) 



where the overline denotes averaging over all particles 
and t w is the "waiting" time, t w < t. 

It is simple to show (see appendix B) that the autocor- 
relation evolves according to the following linear equation 



— A(t w .t) = — A{t w ,t), 



(37) 



where time is again expressed in terms of the collision 
counters t w and r. Equation (|37]) is solved to give 
A(t w ,t) = A(t w ,t w ) exp[-i^(r - t w )], or equivalently 



A(t w ,t) = A (l + t w /Q 1/e 2 {l+t/U) 



(38) 



with A = dT . Therefore, A(t w ,t) is a function of the 
waiting time t w and the observation time i, and not sim- 
ply of their difference, t — t w . This interesting history de- 
pendence or "aging" is another signature of the nonequi- 
librium nature of our system. 

Memory of the initial conditions can be quantified by 
setting t w — 0. Writing A(t) = A(0, t) we arrive at the 
following algebraic decay 



A{t) =A (l + t/U) 



-l/e 



(39) 



In contrast with the temperature which decays with a 
universal law, T(t) ~ i~ 2 , the autocorrelation decays 
with a non-universal law, A(t) ~ t -1 / e . The exponent 
is independent of the dimension. However, it strongly 
depends on the dissipation, and the stronger the dissipa- 
tion, the stronger the memory of the initial conditions. 
This decay exponent is bounded by 2 < f/e < 00. In the 
elastic case, e = 0, a simple exponential decay occurs, and 
in the totally inelastic case, e = 1/2, the autocorrelation 
and the temperature are proportional to each other. 

The autocorrelation function allows calculation of the 
long-time spread in the position of a tagged particle 
A 2 (i) ee (|x(t)-x(0)| 2 ). Using x(t) - x(0) =f*dt'v(t'), 
one can immediately express A 2 (i) via the autocorrela- 
tion function, A 2 = 2 f* dt' /*' dt" A(t" ,t'). Substituting 
(|38| ) into this expression and performing the integration 
yields A 2 (i) = Ci In (1 + t/Q + C 2 [(1 + t/Q 1 - 1 ^ - l] 
with Ci = 2A Q tle/{\ - e) and C 2 = -Ci e/(l - e). 
Asymptotically, the second term is negligible, and the 
spread has a generic logarithmic behavior 

A - Vhii (40) 
reflecting the i _1 decay of the overall velocity scale 

mi- 
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V. STEADY STATES 

Thus far, we have discussed freely cooling systems 
where the energy decreases indefinitely. In typical ex- 
perimental situations, however, the system is supplied 
with energy to balance the energy dissipation |^-g,^3|. 
Theoretically, it is natural to consider white noise forc- 
ing |2^,|^|, i.e., coupling to a thermal heat bath which 
leads to a nonequilibrium steady state. Interestingly, a 
stretched exponential behavior, P[v) oc exp(— v 3 / 2 ), is 
found for the driven inelastic hard sphere gas J22|. 

Specifically, we assume that in addition to changes due 
to collisions, velocities may also change due to an external 



forcing: 



dvj 

TFT 



I heat = £j with j = 1, . . . , d. We use standard 
&(*)&(*')} =2D5 i:j 8{t-t') 



uncorrelated white noise 
with a zero average = 0. The rate equation for the 
temperature is modified by the additional source term 
JjT + AT 3 / 2 = 2D, and the system approaches a steady 
state, Too = (2D/A) 2 / 3 . The relaxation toward this state 
is exponential, (T^ - T\ ~ e ~ const - x *. 

Uncorrelated white noise forcing amounts to diffusion 
in velocity space. Therefore, Eq. (0) is modified as fol- 
lows, -g- — > + Dk 2 . In the steady state, the Fourier 
transform, F^k) = ^S>(y) with y — Dk 2 , obeys 



(i + tf)*(y) = (#(fr)*(w)>. 



(41) 



This equation is solved recursively by employing the cu- 
mulant expansion 



= CX P 



E^(-2/) r 



(42) 



The cumulants K n , defined as 



Foo(k) = cxp 



E 



(43) 



are related to the coefficients ip n , viz. k„ = (2n)lD n ip v 



Writing 1 + y — cxp 
into 



exp 



E n >x(-y) n /n 



we recast Eq. 



(44) 



with the auxiliary variables ip n = V> n (l— £ n — r/ n ). The de- 
sired cumulants i\) n are obtained by evaluating recursively 
the angular integrals of the auxiliary variables, (ip n ), and 
then using the identities ip n = (V>n)/A„. In one dimen- 
sion, (/i™) = 1 and one immediately obtains (ip n ) = n -1 , 
and consequently nip n = [1 - e 2n - (1 - e) 2rl ] 1 [|9). 
In higher dimensions, the quantities (0>n) acquire non- 
trivial dependence on n, e.g., (tpi) = 1, (-02) = ^(V'i), 

and {tp-i) = {ipxip2l — \ The first few values for tp n 
can be then evaluated. In particular, ipi = 1/Ai and 



"02 = ((1 — £, — 7 7) 2 )/(2A 2 A2) , from which one can deter- 
mine explicit expressions: 



1p2 



2e(l-e) 



(45) 



3d 2 



A{d + 2)(1 - e 2 ) - 12(1 - e) 2 (l + e 2 ) ' 



Thus, the steady state distribution is not purely 
Maxwellian. 

To probe velocity correlations or alternatively, devia- 
tions from a factorizing Maxwellian distribution, we con- 
sider the quantity Q, defined in Eq. (p6|). At the steady 
state, it is given by 



Q 



W"(0) 
[*'(0)] 2 



1. 



(46) 



In terms of the first two coefficients of the cumulant ex- 
pansion, Q = 20 2 /0i- Substituting the value of these 
coefficients yields 



6e 2 (l - e) 



(d + 2)(l + e)-3(l-e)(l + e 2 )' 



(47) 



Note that for a fixed spatial dimension, this quantity is 
maximal in the completely inelastic case. For instance, 
Qmax = 2/11 in two dimensions and Q max = 2/15 in 
three dimensions. These values are smaller by an order 
of magnitude or more than the corresponding values in 
the unforced case. Intuitively, one expects that white 
noise forcing randomizes the velocities of the particles. 
Indeed, velocity correlations are much less pronounced 
in this case, as seen in Figure 3. Additionally, veloc- 
ity correlations diminish as the dimension increases. At 
large dimensions, velocity correlations vanish according 
to Q ~ d , indicating that the velocity distribution 
becomes purely Maxwellian, ^(y) — > exp(— y/2), when 
d — ► oo. 



O 



6=2 

d=3 




FIG. 3. The velocity correlation measure Q versus the dis- 
sipation coefficient e. The scaling regime resu lt ( ] 2T[) is shown 
in the top graph, and the steady state result (llTpis shown in 
the bottom graph. 
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VI. GENERALIZATIONS 



VII. CONCLUSIONS 



The above results can be generalized in a number of 
ways. For example, the development of spatial correla- 
tions can be considered by placing particles on a lattice 
and allowing for nearest-neighbor collisions only. In this 
section we briefly mention two straightforward general- 
izations to: (1) energy generating collisions, and (2) dis- 
tribution of restitution coefficients. 

Thus far, we discussed only the physical case of dis- 
sipative collisions, namely e < 0. However, the above 
results in the freely evolving case hold for energy gener- 
ating collisions, i.e., e > as well. Although the typical 
velocity scale diverges, the velocity distribution still fol- 
lows the scaling solution (|Tl]) with algebraic large velocity 
statistics. The corresponding exponent a is still obtained 
from Eq. (|l9|). However, the behavior does change, as fol- 
lows from the analytically tractable d — > oo behavior. In 
contrast with the dissipative case, the second term on 
the right-hand side of Eq. ([l9]) now dominates, and it 
grows exponentially as oc a^~ d ^ 2 . Since the left-hand 
side of Eq. ( p^ ) is of order unity, the constant a must be 
equal to one. On the other hand, the constant a is eval- 
uated using the Stirling formula T(x) ~ {x/e) x to give 
a = (1 - e) 2 (er - d)d d ^ rr -^ a^^-^ . Equating a = 1, 
we arrive at 

o = dv, with {v- \)v~^ = (l-e)~ 2 . (48) 

While the exponent rises linearly with the dimension, it 
exhibit different e-dependence. Numerical solution of v 
shows that this large dimension estimate again yields a 
useful approximation even at moderate dimensions. 

Several recent studies have used a distribution of resti- 
tution coefficients to model driven granular systems, in- 
cluding for example, a one-dimensional gas of rods with 
internal degrees of freedom |5^j5l[], and vertically vi- 
brated layers [ j52f . By tuning the distribution properly, 
one can have a situation were overall, energy is conserved 
as dissipative collisions are balanced by energy generating 
collisions. When the restitution coefficient is drawn from 
the distribution p(e), one simply integrates the collision 
integral in the Boltzmann equation (^]) with respect to 
the measure pie). In one dimension, one can check that 
the scaling solution $(x) = (l + s/x)e^^ still holds, and 
in particular the exponent a — 4 is robust. In general di- 
mension, the exponent a is given by 

I - A ^ = J dep(e)(^- d ^ 2 + 7 /- d V 2 ) (49) 

with the decay rate A = J dep(e)X(e). We conclude 
that algebraic large-velocity statistics extend to situa- 
tions where the dissipation coefficient e is drawn from a 
given distribution. 



We have studied inelastic gases within the framework 
of the Maxwell model, a Boltzmann equation with a uni- 
form collision rate. We have shown that this kinetic the- 
ory is analytically tractable as closed evolution equations 
characterize the Fourier transform and consequently mo- 
ments of the velocity distribution. In the freely evolving 
case, the system approaches a scaling regime, and the 
velocity distribution has an algebraic large velocity tail. 
The corresponding exponent varies continuously with the 
spatial dimension and the degree of dissipation. The de- 
cay exponents can be very large and therefore it may be 
difficult to distinguish a power law from a stretched ex- 
ponential. In the driven case, we have determined the 
cumulants of the velocity distribution. 

The time dependent behavior displays a number of in- 
teresting features. Moments of the velocity distribution 
exhibit multiscaling asymptotic behavior, and knowl- 
edge of the typical velocity is insufficient to character- 
ize all moments. The velocity autocorrelation decays al- 
gebraically with time, and the corresponding exponent 
depends on the restitution coefficient only. 

In contrast with elastic collisions, stochastic inelastic 
collision processes are not effective in mixing particle ve- 
locities. The stronger the inelasticity, the stronger the 
history dependence, i.e., memory of previous behavior. 
Additionally, inelasticity can generate significant corre- 
lations between different velocity components. Such cor- 
relations do develop even in the forced case, where dissi- 
pation is balanced by energy input, and one may expect 
that Maxwellian velocity distribution emerge. 

The Maxwell model is truly mean field in nature with 
all aspects of the collision process being random. While 
it is not surprising that such a theory is solvable, the 
rich structure of the solution is somewhat unexpected. 
For example, the exponent follows from a transcendental 
equation, and can not be obtained from heuristic argu- 
ments or dimensional analysis. Remarkably, even the 
leading asymptotic behavior in the large dimension limit 
remains nontrivial as it involves roots of cubic or tran- 
scendental equations. 

We have explored only the basic characteristics. 
Clearly, one can study higher order velocity correlation 
measures as well as higher order autocorrelations. Fur- 
thermore, the relaxation toward the steady state appears 
analytically tractable. The straightforward analysis is 
cumbersome and it may be useful to expand first the 
solutions in terms of more natural building blocks, e.g., 
orthogonal polynomials. 

We stress that the Maxwell model is exact for stochas- 
tic inelastic collision processes with random collision 
partners and impact angles. It may be applicable in sit- 
uations where an effective stirring mechanism leads to 
perfect mixing. Otherwise, it should be regarded as an 
uncontrolled approximation of the Boltzmann equation. 
Indeed, existing theoretical and numerical studies give 
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little evidence for algebraic tails characterizing inelastic 
gases. The only exception was observed in a system with 
random restitution coefficients drawn from a broad distri- 
bution. In one dimension, both molecular dynamics sim- 
ulation and direct integration of the Boltzmann equation 
for inelastic hard spheres show that the velocity distri- 
bution has a power law tail |5^| . 

In conclusion, our results, combined with previous ki- 
netic theory studies that find exponential, stretched ex- 
ponential, and Gaussian tails, indicate that extremal ve- 
locity characteristics can be sensitive to the details of 
the model, let alone parameters such as the restitution 
coefficient, and the dimension. 

We thank A. Baldassari and M. H. Ernst for fruitful 
correspondence, and H. A. Rose for useful discussions. 
This research was supported by DOE (W-7405-ENG-36), 
NSF(DMR9978902), and ARO (DAAD19-99-1-0173). 
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APPENDIX A: THE A-COEFFICIENTS 

To compute the coefficients A„ = (1 — £ n — rf 1 ) and 
Km = <e^ m )weuseC= l-(l-e 2 )^and7 ? = (l-e)V 
Thus, the following integrals are required 

,„„v = r (l) r ("+l) = 1_3_ 2(n-l) + l 
^ 7 r(|)r(n+|) d2 + d 2(n-l)+d" 

In particular, (/i) = ±, ( M 2 ) = ^^y, ( M 3 ) = , 
so the first few coefficients are 

A 1 =2f(l-e)i, 

A2 = 2(1 _ a l- 2(1 - ()2(l + e ,)_J_, 



+2e(l-e) d (3 + e 2 ) 
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d(d + 2)(d + 4)' 



v(t + At) = 



v prob. 1 — At, 

v — (1 — e)(v — u) • n n prob. At. 



Here u is chosen randomly from all particles and the 
impact direction n is drawn from a uniform distribu- 
tion. The rate of change in the autocorrelation function 
A(t w , t) = v(t w ) ■ v(t) is evaluated as follows 



— A (t w ,t) = lim v(t w ) ■ [v(t + At) - v(t)] /At 
= -(1-e) / duP(u,T) J dn[v(r w )-n] [(v - u) • n] 
= - —j— V ( T «>) ' v ( r ) + —J- J duP(u,T) v(t w ) ■ u 



-A(t w ,t). 



(Bl) 



The angular integration in the second line of Eq. (Bl) 
was performed using the identity 



H(a, b) = J dn (a- n) (b • n) = i (a- b). 



(B2) 



APPENDIX B: THE AUTOCORRELATION 
EVOLUTION EQUATION 

It is useful to work with the collision counter r. In an 
infinitesimal time interval At, the velocity of a particle 
changes from v = v(t) to 



This identity can be deduced by re- writing the integral as 
-ff(a, b) = a • h(b). By symmetry, h(b) = J dnn(h ■ n) 
is a vector along b, say Ab, implying H (a, b) = A (a • b). 
Evaluating the special case if (a, a) = (/i) a 2 we obtain 
Finally, the second term in the third line vanishes, 



v(t w ) ■ u(t) = 0, since the velocity u(t) of the randomly 
chosen collision partner is uncorrelated with v(t w ). 
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